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Abstract 



A spectral parameter power series (SPPS) representation for regular solutions of singular 
I Bessel type Sturm-Liouville equations with complex coefficients is obtained as well as an SPPS 

. representation for the (entire) characteristic function of the corresponding spectral problem on 

a finite interval. It is proved that the set of zeros of the characteristic function coincides with 
the set of all eigenvalues of the Sturm-Liouville problem. Based on the SPPS representation 
' a new mapping property of the transmutation operator for the considered perturbed Bessel 

operator is obtained, and a new numerical method for solving corresponding spectral problems 
is developed. The range of applicability of the method includes complex coefficients, complex 
spectrum and equations in which the spectral parameter stands at a first order linear differential 
operator. On a set of known test problems we show that the developed numerical method based 
on the SPPS representation is highly competitive in comparison to the best available solvers 
such as SLEIGN2, MATSLISE and some other codes and give an example of an exactly solvable 
0^ . test problem admitting complex eigenvalues to which the mentioned solvers are not applicable 

meanwhile the SPPS method delivers excellent numerical results. 

m ■ 
o 

^ ; 1 Introduction 

In the present work the equation 



+ q{x) ] u = X{ri{x)u + ro{x)u) , xe(0, a], (1.1) 



is studied, where / is a real number, I > q is a complex- valued continuous function on (0,a] 
satisfying a growth bound \q{x)\ < Cx" at the origin for some a > —2, ro,i G C [0, a] are complex 
valued functions, and A is a (complex) spectral parameter. Denote L = —-^ + + Qi-i^)- 

Equations of the form (jl.ip appear naturally in many real-world applications after a separation of 
variables and therefore have received considerable attention (see, e.g., [6], [11], [16], [23], [S]). In 
preceding publications equation (jl.ip was studied under more restrictive conditions, typically for 
q and tq being real- valued and ri = 0. The approach developed in this work does not imply such 
restrictions and serves both for qualitative study of solutions and spectral problems, as well as for 
related numerical computation. 

The main component in the developed approach is a spectral parameter power series (SPPS) 
representation for the regular solution of (jl.ip obtained under the condition that the auxiliary 
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equation Luq = possesses a regular solution which does not have zeros on [0, a] except at x = 0. 
The SPPS representation for solutions of nonsingular linear differential equations and its applica- 
tions in corresponding scattering and spectral problems were studied in [9], [10], [15], [T8], [21], 
[22], [23], [25], [26], [19], [27], [32], [36] and some other papers. Here, in Section [2] we obtain an 
analogous result for the perturbed Bessel equation (|1.1|) . The construction and the existence of 
the required particular solution are addressed in Section [3l For example, when q{x) >0,x£ (0, a] 
such nonvanishing on (0, a] solution exists. We give an analytic representation for it together with 
an estimate. Let us emphasize that, in general, uq is allowed to be a complex-valued function, and 
the existence of such uq for a complex valued q is an open problem. 

Under the assumption that uq exists we obtain a dispersion (characteristic) relation for the 
Sturm-Liouville problem for (11. ip on [0, a], a < oo (Theorem 14. ip . Namely, we construct an entire 
function <&(A) in the form of a Taylor series such that the set of its zeros coincides with the set of 
all eigenvalues of the Sturm-Liouville problem. This immediately implies the discreteness of the 
spectrum and offers an efficient method for its computation. 

In practical applications of the SPPS method it is often convenient to consider not only series 
with the centre in the origin A = but also series with the centre at A = Aq where Aq is some 
complex number. In Section [5] we show that this spectral parameter shifting technique is applicable 
to equation ()1.1|) and give necessary details. 

The SPPS representation allows us to obtain a result on mapping properties of the transmuta- 
tion operator corresponding to the operator L, which was studied, e.g., in [Tlj, and [IQ]. We 
show in Section [6] how the transmutation operator acts on certain powers of the independent vari- 
able. In the case of nonsingular Schrodinger operators a result of this kind allowed us to advance 
in the construction of the transmutation operator itself [30j and had applications in constructing 
complete systems of solutions for some partial differential equations [7] , [H] , [20j . 

In Section[7]we explain the numerical implementation of the developed SPPS method for solving 
Sturm-Liouville problems for (II. ip . First, we consider several known test problems, and comparing 
the obtained results with the results obtained by the best available codes, as SLEIGN2, MATSLISE 
and some others, we show that our method is highly competitive and gives better or at least 
comparable results on test problems to which other codes are applicable. Second, we consider 
an example which involves a different from zero ri in (jl.ip and a complex spectrum. Meanwhile 
SLEIGN2 and MATSLISE are not applicable to problems admitting complex eigenvalues, our 
method delivers results which are in excellent agreement with the exact data. 

2 Construction of the bounded solution of a perturbed Bessel 
equation 

Consider a perturbed Bessel operator (also known as a spherical Schrodinger operator) 

(f 1(1 + 1) , , 1 , , , , 

L = --^ + ^^+q{x), l>--,x€{0,a], (2.1) 

where the potential q is (in general) a complex-valued continuous function on (0, a] satisfying the 
growth condition in the origin 

q{x) = 0(x"), X for some a > -2. (2.2) 

Note that we understand the O-notation in the sense of inequality, i.e., there exist a neighborhood 
(0, e] of zero and a constant C > such that \q{x)\ < Cx" for all x G (0,e]. 
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If / 7^ or g L^(0, a], the left endpoint is singular. Despite of that, the equation Lu = 
possesses a solution (j){x) which is bounded at a; = and satisfies the following asymptotics at x = 



(x) ~ X 



l+i 



x^O, 



(p'ix) ~ (/ + X ^ 0, 



(2.3) 
(2.4) 



see, e.g., [241 Lemma 3.2] for the real- valued potential q. We show the explicit construction of 
the solution with this asymptotics at zero for the general case of complex potentials in Section [3] 
meanwhile in Section |3] we show that such solution is unique. 
Together with L consider a linear differential operator 



Ru = vqu + riu' 



(2.5) 



of order at most one, where ro,i G C[0, a] are complex-valued functions, and consider the following 
differential equation involving a spectral parameter A 



Lu = XRu 



or 



u" + 



^^^^J^ + q{x) ) u = X{ri{x)u' + ro{x)u) . (2.6) 



In order to construct a spectral parameter power series representation of a non-singular in zero 
solution of (12. 6p assume that there exists a non- vanishing on (0, a] complex- valued solution uq of 
the equation 

'1(1 + 1) 



+ q{x) uo = 



(2.7) 



satisfying together with its first derivative the asymptotic relations (j2.3p and (j2.4p . Let us define 
the following system of recursive integrals 



1, 



uo{t)R[uo{t)X^''-^\t)] dt, if n is odd, 

^ x("-i)(t) 



(2.8) 



ulit) 



dt, 



if n is even. 



We keep the notation X for consistency with the notations from other publications on the SPPS 
method, see, e.g., [27], [H], [29]. Note that for an odd n we have 



R[uoX(--'^] = {nl+ro){uoX^-'^) 



nu'oX^^-^^ - nuo 



Xin-2) 



Un 



+ roUo 



Uo 



Hence we can write ()2.8|) in a different form, which does not require differentiation of the functions 

''\uo{t)R[uo]{t)X^"-^\t) - ri{t)x'^'^-^\t)) dt, if n is odd. 



x("-^)(t) 



dt. 



The following lemma establishes that all the involved integrals in 
provides some estimates for the functions X^"'\ 



(2.9) 

if n is even. 

are well defined and 
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Lemma 2.1. Let (|2.7p admit a solution uq G C[0,a] nC^(0,a] (in general, complex-valued) which 
does not have other zeros on [0, a] except at x = and satisfies the asymptotic relations (j2.3p and 
(j2.4p . T/ien t/ie system of functions {-^'■"•'}^q is well defined by (j2.8p or ()2.9p and the functions 
satisfy the inequalities 



{2l + 2)r 



+ l)c2n+1^2(/+l)+n 



(2Z + 2; 



n+1 



(2.10) 



where {x)n = ^y'(x)^^ = x{x + 1) ■ . . . ■ (x + n — 1) is the Pochhammer symbol, C = max{l, Ci, C2, C3} 
and the constants Ci, C2 md C3 are such that for any t € (0, a] the following inequalities hold 



\uoit)R[uo]{t)\ < Cit 



21+1 



ul{t) 



< C2t 



-21-2 



\ri{t)\<C^. 



(2.11) 



Remark 2.2. The constants Ci and C2 in Lemma [2 . 1 1 exist due to the fact that uq is a non- vanishing 
function possessing asymptotics (|2.3p and (|2.4p . 

Proof. The proof is by induction. Indeed, for n = we have |X(°)(x)| < 1 and 



X^^\x)\< / \uQ{t)R[uQ\{t)\dt< / Cit 



t.21+1 



dt < 



Cx^ 
21 + 2' 



Assuming that the statement is true for some n = k, for n = k + 1 we obtain 



|l(2(fc+l))(^)| < 



X(2fc+l)(t) 



ulit) 



dt < 



C2 {k + i)C^k+i^2{i+i)+k 



l2l+2 



{21 + 2)k+i 



dt < 



(j2k+2^k+l 
{2l + 2)k+l 



and 



^(2(fc+iHi)(^)| < r\uo{t)R[uomx^^'+^'>{t)\dt+ r\n{t)X^^''+'\t)\dt 
Jo Jo 



< 



< 



Cit 



21+1 



c 



2k+2^k+l 



{21 + 2)fc+i 

(j2k+3^2l+2+k+l 



dt + 







+ 



^ {k + l)(72fc+1^2(/+l)+fc 
Cq • ; dt 

(2/ + 2)fc+i 

C^k+2 . ^ -^)^2/+2+fc+l 



< 



(2Z + 2 + A: + 1) • (2/ + 2)^+1 (2/ + 2 + A: + 1) • (2/ + 2)k+i 
{k + 2)c2'=+32;2a+i)+fc+i 



{21 + 2)fc+2 

Note that the exponents of the powers of t in all the involved integrands are non-negative, hence 
all the recursive integrals are well defined. □ 

In the particular case when ri = 0, i.e., the right hand side of (j2.6p does not depend on the 
derivative of u, the estimates of Lemma l2. II can be improved, and the following statement is valid. 

Lemma 2.3. Under the conditions of Lemma \2.1\ assume additionally that ri = 0. Then the 
functions X*^"^ defined by (12. 8p or (12. 9p satisfy the inequalities 



|^^^")(x)| < 



2ra^2n 



C72"x 



22"n!(/ + 3/2)„' 



|x(2"+i)(x)| < 



(j2n+1^2n+l+2(l+l) 



22'"+in!(/ + 3/2) 



n+1 



(2.12) 



where {x)n is the Pochhammer symbol and C = max{Ci,C2}, where the constants Ci and C2 are 
such that for any t £ (0, a] the following inequalities hold 



\rQ{t)ul{t)\ < Cit 



21+2 



1 



ul{t) 



< Cot'^^'^. 



(2.13) 
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The proof is analogous to that of Lemma 12.11 

The following theorem presents the spectral parameter power series (SPPS) representation of a 
bounded solution of equation (j2.6p . 



Theorem 2.4. Let (|2.7|) admit a solution uq G C[0, a]nC^(0, a] (in general, complex-valued) which 
does not have other zeros on [0, a] except at x = and satisfies the asymptotic relations (j2.3p and 
(j2.4|) . Then for any A S C the function 



oo 
k=0 



uoX^A'^X^^fc) (2.14) 



is a solution of (j2.6p belonging to C[0,a] n C^(0, a] and the series converges uniformly on [0,a] 
The first derivative of u is given by 



/ ^0 
u = —u 



Lyx^x'^^k-i)^ (2.15) 



no no 

and the series for the first and the second derivatives converge uniformly on an arbitrary compact 
K C (0,a]. 

Proof. Formally differentiating the series (j2.14p twice with the aid of (j2.9p we obtain that u' should 
be given by (j2.15p and u" (after simplification) by ^u— Xr^u— Xriu' . By Lemma [2. II all the involved 
series converge uniformly on an arbitrary compact K C (0, a], hence the formal derivatives coincide 
with the usual ones. Since by (12. Sp and (12. 7j) 



Uo V ^ / 



+ q{x) I u — XR[u], 



u is indeed a solution of equation (j2.6p . The relations (j2.3p and (j2.4p follow from the corresponding 
asymptotics of uq because by Lemma 12.11 we have X(2fc)(x) = o(l) and ^^(^'^-i) = o(xO for 
k>l. ° □ 

Remark 2.5. For a regular Sturm-Liouville problem the existence and the construction of the 
required non-vanishing solution uq present no difficulty since the equation possesses two linearly 
independent real-valued solutions ui and U2 whose zeros alternate and one may choose uq = ui+iu2 
as such solution. For the singular equation under consideration there is only one solution satisfying 
the asymptotic conditions (|2.3p and (j2.4p , see Theorem 14.11 Corollary 13.31 establishes that such 
non- vanishing solution uq exists in the case when q{x) > 0, x G (0, a], and in Remark 15.21 we 
show that a modified SPPS representation is always possible in the case when q is real valued and 
bounded from below and ri is real valued. 

Remark 2.6. For SPPS representations for solutions of nonsingular Sturm-Liouville equations we 
refer to |25], [26] and [27] . They have been applied in a number of papers to different scattering 
and spectral problems (see references in the Introduction). For the perturbed Bessel equation, in 
the case of a real valued potential g, ri = and tq = 1 an SPPS representation was obtained and 
used in [24] but without formulas for constructing or estimating the coefficients X^'^^\ 

For practical applications the partial sums of the series ()2.14p are of the main interest. Based 
on Lemmas 12.11 and 12.31 the following corollary provides estimates for the difference between the 
exact solution and the approximate one defined as a partial sum of the series (j2.14p . The difference 
is estimated in terms of the remainders of Taylor series of two special functions. 



5 



Corollary 2.7. Under the conditions of Theorem \2.4\ consider um = uo'^k=o ^^■^^'^^'^ ' '^ote that 
for N = the right-hand side is equal to uq. Then 



\u{x) - un{x)\ < max \uo{t)\ ■ ^ ' ' < max \uo{t)\ ■ 1^1^ - ^ 



N 

I ( \ Al-r 1"- 

2.16) 



kl 

k=0 

where the constant C is defined in Lemma \2.1\ 

Moreover, in the particular case ri = the following estimate holds 

|n(x)-n^(x)|< maxMOl- ^ 22fc^!(/ + 3/2), 

(|A|i/2Cx)'+i/2''+i/2^l^l A.22feA;!(/ + 3/2)fe 



< max |'Uo(t)| 

te[0,a;] 



where r(/ + 3/2) is i/ie Gamma function, Ii^i/2{x) is the modified Bessel function of the first kind 
and the constant C is defined in Lemma \2.3i . 



For the proof one should simply compare the majorizing terms corresponding to the even in- 
dices in ()2.10p and (|2.12p with the Taylor expansions for the exponential and the Bessel functions 
appearing in (j2.16p and (|2.17p . 

Example 2.8. Consider the Bessel equation 

d^ Kl + 1) 
— -r-^u H K — u = \u, 

with / > —1/2 and x G (0, a]. The regular solution of this equation satisfying ()2.3p and ()2.4p is 
given by the formula 

ui{x, A) = r(/ + 3/2)2^+1/2 A-"^ . ^J^^,/2(^^), (2.18) 

where Jl-^-l/2 is the Bessel function of the first kind. This solution may be represented as a power 
series in terms of the parameter A, 

00 { -\\k 2k 

see 



the equation —u" + u = 0, satisfying (j2.3p and (j2.4p . It is easy to verify that choosing such 



In order to apply Theorem 12.41 consider uq{x) = x'^^ as the non- vanishing on (0,a] solution of 
; equation —u" + ^ 
uq we obtain from (12. 9p 

^ ( -\\nrJ2n 



4"n!(/ + 3/2)„' 
i.e., exactly the coefficients from (|2.19p . 

3 Construction of the particular solution uq 

In this section we explain how to construct a particular solution of equation ()2.7p satisfying asymp- 
totics (|2.3p . (|2.4p and present some sufficient conditions for this solution to be non- vanishing for 
X > 0. 
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In order to construct an SPPS representation for the particular solution we rewrite equation 
(j2.7p in the form 

y" - ^^^y = q{x)y. (3.1) 



The equation 



possesses two solutions 



y'o - ^-^^yo = (3-2) 



yi{x) = X ' and ^2(2;) = 2;'^^. 

For / > —1/2 the second solution is regular and satisfies (j2.3|) and (|2.4p . Since the potential q may 
be singular in the origin, we cannot apply Theorem 12.41 directly. However we may construct the 
system of recursive integrals in the same way as in (j2.9p and only have to show the convergence of 
the integrals and obtain some estimates justifying the SPPS representation. 
Consider the following system of recursive integrals 

= 1, 



y("-i)(t)t2('+i)g(t)dt, for odd n, (3 3) 

/ y("-^)(t)t-2('+i)dt, for even n. 
'0 

Note that since the potential q G C(0, o] satisfies condition (j2.2p for some a > —2 there exists 
a constant C > such that 

\q{x)\ < Cx" for ah x G (0, a]. (3.4) 

Lemma 3.1. Suppose that the complex-valued potential q G (7(0, a] satisfies inequality ()3.4p for 
some C > and a > —2. Then the functions y(") are well defined by ()3.3p and the following 
estimates hold 

l^'"'<"*l ^ (2 + o)-„!(|±l + l)„' 

_ fin 2l+l+n(2+a) 

l^'^"-"'-'l ^2.o,^..-„,(|S.l)„ ' -<»-l' P-^' 

where (x)„ is the Pochhammer symbol. 

Proof. The proof can be performed by induction, similarly to the proof of Lemma 12.11 The only 
difference is that it is possible that an exponent of the power of t under the integral sign be negative, 
however it is always strictly greater than —1. Hence all the involved integrals exist. □ 



Proposition 3.2. Suppose that the complex-valued potential q G C(0,a] satisfies inequality ([37 
for some C > and a > —2. Then the function 



00 



no(x) = x'+^^y(2^')(2;), (3.7) 



k=0 



where the functions Y^'^^^ are defined by (|3.3p . is a particular solution of equation (j3.ip on (0, a] 
satisfying asymptotics (12. 3p . (12. 4p . Moreover, uq satisfies the following estimate for any x G (0, a] 

/2; + l \ 21+1 21+1 ^ ( 2VCx2+«\ 

\uo{x)\ < T + 1 (2 + a) — "1 . (3.8) 

\ 2 + a / 2+a \ I -\- a j 
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Proof. Due to ()3.3p the first and the second derivatives of uq are given by the expressions 



[, = (/ + J2 Y^^'^ + E ^^'""'^ (3-9) 



k=0 k=l 



and 



2fe-l) 



k=0 k=l 
k=l k=l 

The uniform convergence of ah the involved series on an arbitrary compact K C (0, a] and hence 
the possibihty of termwise differentiation, fohows from Lemma |3.1[ The asymptotic relations (j2.3p 
and (|2.4p for the function uq follow from the estimates y^^^) = o(l) and y(2fc-i) = o(x^'+-'^), x — )• 
for A: > 1, see (j3.5p and (|3.6p . Inequality p.Sp follows from the series representation of the modified 
Bessel functions of the first kind, see, e.g., [ij. □ 

The following corollary provides a sufficient condition for the particular solution constructed in 
Proposition 13.21 to be non- vanishing. 

Corollary 3.3. Under the conditions of Proposition \3.2\ assume additionally that q{x) > 0, x G 
(0, a]. Then Uo{x) > x'"*"-^ for any x S (0, o]. 

Proof. We obtain from (|3.3p and the condition q(x) > that Y^^\x) > 0, n > 1. Hence, no(x) = 



4 Spectral problems 

The classical formulation (see, e.g., [24], f42]) of a spectral problem for a Sturm-Liouville equation 
of the form Lu = Xru with L from (j2.ip consists in finding the values of the spectral parameter for 
which there exists a solution u{x] A) continuous at x = and satisfying the following conditions. 
When Z > 1/2, 

m(0;A)=0 (4.1) 

and 

pu{a; A) + 7n'(a; A) = (4.2) 

for some /3, 7 G C such that + I7I / 0. 

When I £ [—1/2,1/2) and since the second linearly independent solution is also square-integrable, 
an additional boundary condition 

lim x' ((/ + l)u(x; A) - xu'(x; A)) = 0, (4.3) 

is imposed (see, e.g., |24]). 

Despite the fact that the right-hand side ()2.5p of the spectral equation ()2.6p may contain the 
derivative of the unknown function and hence not fit into the basic Sturm-Liouville scheme, we 
consider the same spectral problem (|4.ip - (j4.3p for equation ()2.6p . The following statement gives 
us a characteristic equation of the spectral problem under the condition that an appropriate non- 
vanishing solution Uq of (j2.7p exists. 
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Theorem 4.1. Let (j2.7p admit a solution uq G C[0, a] n C^(0,a] ('in general, complex-valued) 
which does not have other zeros on [0, a] except at x = and satisfies the asymptotic relations ()2.3p 
and (|2.4p . T/ien the eigenvalues of the problem (|2.6p . ()4.ip - (j4.3p coincide with zeros of the entire 
function 



oo 



where the functions X^''^ are defined by (|2.8p or (I2.9p . 

Proof. Under the condition of the theorem, the function u defined by (|2.14p is a solution of (j2.6p 
satisfying the boundary conditions ()4.ip and ()4.3p . The boundary condition (|4.2p for the function 
u coincides with <^*(A) = 0, and ^{X) is an entire function by Theorem 12.41 It is left to show that 
there are no more eigenvalues. For that it is sufficient to show that the second linearly independent 
solution U2 of (j2.6p does not satisfy either (j4.ip or (j4.3p . We rewrite (j2.6p in the form 

- u" - Arm' + C-^!-^ +q- Xro) u = (4.5) 



and introduce the function 

p{x) = e-^-^o'^^^^^'^K 

Note that 

lim p{x) = 1. (4-6) 

x—^0 

Since the solution u defined by (|2.14p satisfies the asymptotic condition (j2.3p . there exists a constant 
b = 6(A) > such that u{x) 7^ for all x S (0, b]. Then the second linearly independent solution is 
given by the Liouville formula [171 Chap. XI], [37J 

U2(x) = —u(x) [ dt, X £ (0,b]. 

Jx U''{t) 

Assume first that I > —1/2. It follows from (j4.6p . asymptotics (|2.3p . ()2.4p and L'Hospital's rule 
that ^ 

^2(3;)--^^, x^O. (4.7) 

Hence for / > the second solution U2 does not satisfy the boundary condition (j4.ip . Let now 
— 1/2 < / < 0. Recall that by Abel's identity [171 Chap. XI] the Wronskian of u and U2 has the 
form 

W = UU2 — UU2 = p. (4-8) 
Hence using (fO]) . ([23]), (fOI) and (gZl) we obtain 

~ x^O, (4.9) 

and thus observe that U2 cannot satisfy the boundary condition ()4.3p . 
For / = —1/2 similar reasoning shows that 

ii2(x) ~ -a/xIux, X ^ 0. (4.10) 

We substitute the expression for u'2 obtained from (j4.8p into ()4.3p and obtain 

-1/2 /"l /\ n-23;n' 
X / -n2 - XM2 = -Va;- H „ r- ■ 
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For the first term in tliis expression, from (j4.6p and ()2.3p we have 



Hm ^/x?^ ^ 



x^o u{x) ' 

hence to prove that U2 does not satisfy the boundary condition (j4.3p it is sufficient to show that 
the second term is o(l) as x — >• 0. Due to the asymptotic relations (I2.3p and (I4.10p it is sufficient 
to show that u — 2xu' = 0(x^/^+^) for some e > 0. Taking into account (|2.15p we have 

u - 2xu' = — (no - 2xu') + \^ X^''^-^\ 

As can be seen from (|2.10p . ^ X^fcLi A'^X^^^"^^ = o(x), hence only the first term is relevant. Since 
u ~ no ~ \/x as X — )• 0, it is sufficient to prove that no — 2xnQ = 0(x^/^"'"'^). Similarly to (I4.10p we 
obtain that the general solution of equation ()2.7p can be represented as ciui + C2U2, where ui is 
given by ()3.7p and satisfies the asymptotic condition (12. 3p and U2 satisfies the asymptotic condition 
U2 ~ ^/x\nx, X — )• 0. Since no ~ \/x by the statement of the theorem, it is necessarily of the form 
(I377l> . Now using ([32D, and ([3S]) we obtain 

00 „ 00 00 

no - 2xn'o = V^Y.^^''^ " ^E^^'"^ " = 

k=0 ^ k=0 k=l 

where a participates in (j2.2p . and this finishes the proof for / = — 1/2. □ 



5 Spectral shift technique 

The SPPS representation given in Theorem 12.41 is based on a particular solution of equation ()2.6p 
for A = 0. In [27J it was mentioned that for a classic Sturm-Liouville equation it is also possible to 
construct the SPPS representation of a general solution starting from a non-vanishing particular 
solution for some A = Ao. Such procedure is called spectral shift and has already proven its 
usefulness for numerical applications [27], [18] . 

We show that a spectral shift technique may be applied to equation (j2.6p . Let Aq be a fixed 
complex number. We rewrite ()2.6p in the form 



u" - Xoriu' + (^-^—^-^ + Q - ^0^0^ u = \{riu + ron), (5.1) 



where A := A — Aq. Suppose that no is a solution of the equation 



Lon := -u" - Xonu' + ( iii±il + ^ _ Aoro ) n = (5.2) 



such that no does not vanish on (0, a]. Note that no is a particular solution of (|2.6p for A = Ao- 
Then the operator Lq admits the following Polya factorization |35] 



1 d 2 d u 

Lqu = l-P«oi > (5-3) 

puQ ax ax no 

where 

p(x) = e^o^o"^i(*)'^^ (5.4) 



10 



Based on the factorization (j5.3p we introduce the following system of recursive integrals 

= 1, = 0, 

p(t)uo(t)i2[Mo(t)^^"~^Hi)] dt, if is odd, 

/ — , , o/ s rf^i if is even, 

/o P(t)^^o(0 

Similarly to ()2.9p the recurrent relation (I5.5p can be rewritten as 



(5.5) 



(lj(t)no(t)ii[no](t)Z("-i)(t) - ri(t)Z("^2)(^)^ -f ^ ^^^^ 
zW(x)= <j-^o r.Zin^Du. (5.6) 

— / — , , o , , dt, if n is even. 

I io PitWoit) 

Theorem 5.1. Let (j5.2p admit a solution uq £ C[0, a]nC^(0, a] ('m general, complex-valued) which 
does not have other zeros on [0, a] except at x = and satisfies the asymptotic relations (12. 3p and 
(j2.4p . T/ien /or any A G C i/ie function 



n = noJ](A-Ao)'^z(2fc) (5.7) 

A;=0 

is a solution of (|2.6p belonging to C[0,a] nC^(0,a], and i/ie series (|5.7p converges uniformly on 
[0, a] . T/ie series for the first and the second derivatives converge uniformly on an arbitrary compact 
K C (0, a] and the first derivative has the form 

^' = ^ _ ^ y (A - X^fz^^^-^^ = u'^ + y (A - Ao)' ( ^(2^=) - f ) . (5.8) 

Proof. Note that the function p given by (j5.4p is continuous and non- vanishing on [0, a]. Similarly 
to the proof of Lemma [2. II we see that the functions Z^") satisfy estimates (|2.1Up with the constant 

C = max{l,Ci,C2,C3}, (5.9) 

where 

\p{t)u,{t)R[uo\{t)\ t^'+^ ^ I 

Ci = sup -757— j C2 = max , ^ ^ 9, s 1 1 ^3 = max ri(t) . 

te(o,a] te[o,a] |p(t)n^(t)| te[o,a] 

Formal application of the operator Lq to (|5.7p with the use of Polya factorization (|5.3p and formula 
(j5.5p shows that the function n is a solution of equation (|5.ip and hence of (|2.6p . Estimates for 
^(") justify the possibility of differentiation of the involved series and show that the function u 
satisfies (EISl) and (El). □ 



Remark 5.2. Suppose that the functions q and ri are real-valued and that the potential q is bounded 
from below, i.e., there exists a constant go £ such that 

q{x) > qo for all x S (0, a]. 

Consider equation (j5.2p for Xq = qo- The particular solution of this equation can be constructed 
similarly to Section [3] using the generalization of formulas (jS.Sp according to the factorization (j5.3p , 
cf., (12. 8p and (j5.5p . Since the functions p and q — qo are non-negative, similarly to Corollarv 13.31 
we deduce that the particular solution no in such case does not have other zeros on [0, a] except at 
X = 0. Hence it is possible to construct the SPPS representation of the bounded solution for any 
equation (|2.6p having a real-valued ri and a real-valued bounded from below q. 
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Remark 5.3. In the case when ri = 0, we do not need to introduce the new system of functions 
^("). The system of functions X^"^ given by (j2.9p may be used directly in representation (j5.7p . 

Remark 5.4. For the difference max^g[o,a] \u{x) — un{x)\, where un = uq XlfeLo ^^Z^'^''^ the estimate 
()2.16p holds, where C is given by ()5.9p . 



6 Transmutation operators for perturbed Bessel operators 

We recall a general definition of a transmutation operator from [29] which is a modification of the 
definition given by Levitan [33j. Let ii^ be a linear topological space and Ei its linear subspace (not 
necessarily closed). Let A and B be linear operators: Ei — )• E. 

Definition 6.1. A linear invertible operator T defined on the whole E such that Ei is invariant 
under the action of T is called a transmutation operator for the pair of operators A and B if it 
fulfills the following two conditions. 

1. Both the operator T and its inverse are continuous in E; 

2. The following operator equality is valid 

AT = TB 

or which is the same 

A = TBT-\ 

Very often in literature the transmutation operators (the term coined by Delsarte and Lions 
|14j ) are called transformation operators. 

In [8] for the case of the transmutation operator T corresponding to the pair of operators 
A = —-^ + q{x) and B = — ^ a mapping property of T was found. It establishes what is the 
result of action of T on the powers of the independent variable. This mapping property already 
found many applications in the proofs of completeness of infinite systems of solutions of some linear 
partial differential equations and in solving initial and spectral problems, see [7], [8], [20j, [28j, 
|30j . |31j . We present an analogue of the aforementioned mapping property for the transmutation 
operator for the pair of operators A = —-^ + q{x) + -^^r^ and B = —-^ + "^^^j constructed in 
[39], [40]. 

We recall some results from [39], [IQ]. Under certain additional conditions on the potential g, 
discussed below, a bounded solution of the equation 

-y"+{q{x)+^-^^)y = \y (6.1) 

can be represented in the form 

rx 

y{x, A) = ji+i/2{x, A) + / K{x, t)ji^i/2{t, A) dt, (6.2) 

Jo 



where ji+i/2{x-, X) = \/ xy/\Jij^ii2{x\/\) is a solution of the equation 



y" + ^-^^y = Ay (6.3) 
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and Ji-i-i/2 is the Bessel function of the first kind. The integral kernel K is the solution of the 
partial differential equation 



d'^K(x,t) 1(1 + 1) N r ^r.^ ^ d'^K(x,t) 1(1 + 1),^, , 
- ^^^(^' - l(^)K{x, t) = ' - ±-^K{x, t) 

satisfying the boundary conditions 

dR{x,x) ^ 1 lim K{x,t) ■ = 0. 



Moreover, the integral kernel K satisfies 



sup 

0<x< 



/ \K{x,t)\'^dt <oo. (6.4) 

a Jo 



We denote the operator defined by (|6.2p as T. The existence of such operator was established in 
|40j for the case when q is a continuous function on [0, a] and in |39j for the case when / is an integer 
and g is a real-valued function satisfying the condition Jq t"^\qit)\ dt < oo for some < m < 1/2. 
It is mentioned in [38] that the results of [2] allow one to extend the existence of the operator T 
onto arbitrary real values of the parameter I. 

The solution y{x, A) in (j6.2p differs from the solution u{x, A) satisfying the asymptotic condition 
(j2.3p by the factor 2'+^/^2^^i+3/2) ' ^H]' '^^^ series expansion of the function ji+i/2{x, A) is 

j,,y,ix, A) = c..-+'+S where c. = c.(A) = r(fc + i)r(fc + / + 3/2)2^^-+'+V2 ' (6-5) 



k=0 



Similarly to [8, Theorem 7] we substitute the functions y{x,X) and J/+i/2(2;,A) in (j6.2p by their 
series expansions (j2.14p and (|6.5p and obtain 

\(l+l)/2 °° °o px / ^ \ 

^ ' J- I' + 2J fc=o fc=0 ^ k=0 ^ 

= 5^ Cfc (^x^'^+'+i + j (6.6) 



A:=0 

^ r(A: + l)r(A; + / + |)22'=+'+V2 ^ ^ 



The function i^(x, •) is square-integrable on [0, x], see (j6.4p . and the function ii+i/2(-, A) is the limit 
of the uniformly convergent series (16. Sh . hence the possibility to change the order of summation 
and integration in ()6.6p follows from the continuity of the scalar product in ^^[0, a]. 
Since the equality in (16. 6p holds for all x and A, we finally obtain that 

T[3.2fe+/+i] = (_i)fc22'^A:! f / + ^) no(x)x(2^) (x). (6.7) 

V 2 / fc 



7 Numerical implementation and examples 

Based on the results of the previous sections we can formulate a numerical method for solving 
spectral problems for perturbed Bessel equations. 



13 



1 . Find a particular solution uq of equation ()2.7p satisfying the asymptotic conditions (j2.3p and 
()2.4p . Note that due to the proof of Theorem 14.11 it is sufficient to check that the solution 
satisfies only the asymptotic condition (j2.3p . If an analytic expression for the particular 
solution is unknown, one can use a numerical approximation suggested by Proposition 13. 2[ 

2. Check that the particular solution uq obtained in step 1 is non- vanishing for x G (0, a]. Under 
the conditions of Corollary 13.31 it is always the case. If the particular solution has zeros on 
(0,a], the spectral shift technique described in Section [5] can help, either directly as described 
in Remark 15.21 or by finding a suitable value of Ao, complex in general. 

3. Use partial sums of the series ()2.14p and (|2.15p (or, correspondingly, (|5.7p and (|5.8p ) to obtain 
a polynomial 

$jv(A) = (/?no(a) + 7n[,(a)) +^ A'^ (/3no(a) + ju'^ia)) X^^^''\a) - -^X^^''~'\a) (7.1) 

fc=i ^ '"oiaj / 

approximating the characteristic function (|4.4p . 

4. Find zeros of the calculated polynomial <I>7v(A). 

5. To improve the accuracy of the eigenvalues located farther from the point Aq, the spectral 
parameter corresponding to the current particular solution, perform one or several steps of 
the spectral shift technique. 



Since the characteristic function ^{X) is analytic, the Rouche theorem from complex analysis, 
see, e.g., ^2\, provides the stability of the numerical procedure. Indeed, let F be an arbitrary simple 
closed contour on which $ does not vanish. Then if the absolute error of approximation — <I>| 
is less than min|<I>7v| on F, due to Rouche's theorem the functions and (^> — ^jy) + = ^ 
possess the same number of zeros inside F. Hence the procedure described above does not produce 
additional (or on the contrary less) zeros whenever approximates well enough the function 
We illustrate this point below, in Example 17. 7[ 

Before considering numerical examples let us explain how the numerical implementation of the 
SPPS method was done in this work. All the calculations were performed with the help of Matlab 
2010 in the double precision machine arithmetics. The formal powers X^"'^ were calculated using 
the Newton-Cottes 6 point integration formula of 7-th order, see, e.g., [13], modified to perform 
indefinite integration. We choose M equally spaced points covering the segment of interest and 
apply the integration formula to overlapping groups of six points. It is worth mentioning that 
for large values of the parameter I a special care should be taken near the point 0, because even 
small errors in the values of after the division by Uq ~ 2;2(/+i) j^g^^^ large errors in the 

computation of X^^"-) on the whole interval [0,a]. To overcome this difficulty we change the values 
of in several points near zero to their asymptotic values. This strategy leads to a good 

accuracy. The computation of the first 100-200 formal powers proved to be a completely feasible 
task, and even for M being as large as several millions the computation time of the whole set of 
formal powers is within seconds. In the presented numerical results we specify two parameters: N 
is the degree of the polynomial <^Ar in (17. ip . i.e., the number of the calculated formal powers is 2N, 
and M is the number of points taken on the considered segment for the calculation of integrals. 

Let us stress that the formal powers do not depend on the spectral parameter and once calculated 
can be used for computing the solution and/or the characteristic function $(A) for thousands of 
different values of the spectral parameter A without any additional significant computation cost. 
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n 


A„ (SPPS) 


A„ (SLEIGN2, (H) 


Ay, (Exact) 


1 


12.1871394680951 


12.187139459 


12.1871394680951 


2 


44.257559403500 


44.257558912 


44.257559403502 


3 


96.07160483898 


96.071604502 


96.07160483884 


4 


167.62571241787 


167.625711908 


167.62571242058 


5 


258.91930037169 


258.919292439 


258.91930035744 


6 


369 95220905860 


369 952209262 


369 95220926235 


7 


500.72440133519 


500.724370471 


500.72438147579 


8 


651.23517308180 


651.235865279 


651.23579210254 


9 


821.50506498326 


821.486428982 


821.48642898238 


10 


988.97560762340 


1011.476285608 


1011.47628560802 



Table 1: The eigenvalues from Example 17.11 for c = calculated with = 40 and M = 5 • 10^. 

There exist several computer codes for the solution of singular Sturm-Liouville problems. We 
compare our results with the results produced by SLEIGN2 [5] and MATSLISE [33] . Both packages 
can reliably solve a variety of spectral problems for regular and singular Sturm-Liouville problems 
and are considered as basic comparison tools, the second package to our best knowledge is one 
of the most accurate. Both packages work with double precision machine arithmetics and were 
used with the parameters for the highest possible accuracy goals. It should be mentioned that 
the applicability of both packages is restricted to the self-adjoint situation and they do not permit 
neither complex coefficients nor a first order differential operator at a spectral parameter. Thus, to 
illustrate the performance of the proposed method in a situation when a comparison to the existing 
software is impossible, in Example 17.71 we consider a problem for which an exact characteristic 
equation can be written down and compare the obtained numerical results to those calculated from 
the exact equation. 

7.1 Real spectrum 

Example 7.1. Our first numerical example is the Bessel equation, which we already touched on in 
Example 12.81 Consider the following spectral problem [U Example 2] 

-y" + ^y = Ay, o<x<i, 

y(l,A) = 0. 

In this and in all other considered examples the solution y{x,X) should also satisfy the boundary 
conditions (14. ip and (I4.3p at 0. As follows from ()2.18p . the eigenvalues of the problem ()7.2p coincide 

with zeros of the Bessel function Ju{s), where s = \/A and v = \J c + ^. 

We compare the results delivered by the SPPS method to those from [H Example 2] for a 
particular case c = For the SPPS representation ()2.14p the exact particular solution uq = x^^'^ 
was used. The results are presented in Table [1] together with the exact eigenvalues calculated as 
squares of zeros of ^3/4(5) and computed using the Matlab routine besselzero .m by Greg von 
Winckel. 

As can be seen from Table [H the accuracy of the higher eigenvalues calculated by the SPPS 
method decreases. To improve the accuracy the spectral shift technique described in Section [S] is 
applied. We choose that the values of the parameter Aq change along a line in the complex plane 

(n) 

and are given by Aq = 50n -|- 2ni, n = 1 . . . 4000. On each step we compute the solution in the 
next point Ag"^^^ and use this solution as a particular solution for the next step. Since the spectral 
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n \n (exact/MATSLISE) An (SPPS) Ap used (SLEIGN2) 



1 


12, 


,1871394680951 


12, 


,1871394680951 


l+0.5i 


12 


.187139459 


2 


44, 


,257559403502 


44, 


.257559403500 


50+2i 


44, 


.257558912 


3 


96, 


.071604838843 


96, 


.071604838834 


100+4i 


96, 


.071604502 


4 


167, 


,62571242058 


167, 


.62571242056 


150+6i 


167, 


.625711908 


5 


258, 


,91930035744 


258, 


.91930035742 


250+lOi 


258, 


.919292439 


6 


369, 


,95220926235 


369, 


.95220926232 


350+14i 


369, 


.952209262 


7 


500, 


,72438147579 


500, 


.72438147575 


500+20i 


500, 


.724370471 


8 


651, 


,23579210254 


651, 


.23579210250 


650+26i 


651, 


.235865279 


9 


821, 


,48642898238 


821, 


.48642898235 


800+32i 


821, 


.486428982 


10 


1011, 


,47628560802 


1011, 


.47628560801 


1000+40i 


1011, 


.476285608 


30 


8956, 


,5077203636 


8956, 


.5077203638 


8950+358i 


8956, 


.507721371 


50 


24797, 


,222775294 


24797, 


.222775296 


24800+992i 


24796, 


.866878938 


75 


55701, 


,421553437 


55701, 


.421553167 


50000+2000i 


55701, 


.328815911 


100 


98942, 


,625835 


98942. 


,625812 


100000+4000i 


98943, 


.253862034 



Table 2: The eigenvalues from Example 17.11 for c = calculated with the use of the spectral shift 
technique with TV = 40 and M = 5 • 10^. 



problem has a purely real spectrum, among the zeros of the approximating polynomial the ones 
with a small imaginary part are chosen and those which are closer to the number Re Aq than to 

(k) 

any other of the numbers ReAg , k ^ n are stored. The obtained results are presented in Table 
[2] together with the used values of Aq for the spectral shift and the exact eigenvalues. In the first 
column of the table we combine the results produced by the MATSLISE package with the exact 
eigenvalues because in this case all the presented digits coincide. 

As can be seen from Table [21 the spectral shift technique allows us to significantly improve the 
obtained results for the higher eigenvalues. From now on we present only the results obtained with 
the help of the spectral shift technique and specify the implemented rules for choosing spectral 
shifts. 

Example 7.2. The second example is the Boyd equation. Consider the following spectral problem 
^ Example 3] 

i-y"-ly = Xy, 0<x<l 
\y{l,X)=0. 

This equation fits into the general scheme if one takes / = 0. We take the function uq = ^Ji{2^/x) 
as a particular solution of (|2.7p and calculate eigenvalues using the SPPS method with = 40 
and M = 50000. For problem ()7.3p the characteristic equation is known and is given by 

kMk,i/2{2iV\) = 0, (7.4) 

where is the Whittaker function and k = k{\) = (2i\/A)~^, see [U Example 3]. In Table 

[3] we present the obtained results together with the exact eigenvalues calculated from ()7.4p and 
the results obtained using SLEIGN2 and MATSLISE software. Note that some eigenvalues caused 
problems for the MATSLISE package, despite the excellent accuracy of the rest of the results. To 
simplify the reading, in the presented numbers we truncated some digits which do not coincide 
with the correct ones. Relative errors of the first 50 eigenvalues obtained by SPPS, MATSLISE 
and SLEIGN2 are presented on Figure [TJ We emphasize a relative stability of the numerical results 
delivered by the SPPS method in comparison to those computed by MATSLISE and SLEIGN2. We 
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n 


A„ (SPPS) 


Xn (exact) 


A„ (SLEIGN2) 


Xn (MATSLISE) 


1 


7. 


,3739850151752 


7, 


.3739850151751 


7.37398499 


7.3739850151751 


2 


36, 


,3360195952325 


36, 


.3360195952318 


36.33601959522 


36.3360195952318 


3 


85, 


,292582094149 


85, 


.292582094137 


85.29258240 


85.292582094137 


4 


154. 


,098623739770 


154. 


,098623739767 


154.098623739742 


154.098623739767 


5 


242. 


,705559362903 


242. 


,705559362911 


242.705559362924 


242.705559362911 


6 


351. 


,091167129407 


351. 


,091167129418 


351.09116712937 


351.091167129418 


8 


627. 


,155044324547 


627. 


,155044324564 


627.155033 


627.155044324564 


10 


982. 


,239093680177 


982. 


,239093680188 


982.239069 


982.239093680188 


13 


1662. 


,98063088581 


1662. 


,98063088578 


1662.9806308832 


1662.76 


20 


3942. 


,42966385096 


3942. 


,42966385102 


3942.18 


3942.42966385102 


28 


7732. 


,02180519217 


7732. 


,02180519214 


7732.021805177 


7729.4 


30 


8876. 


,82700072940 


8876. 


,82700072941 


8876.8270009 


8876.82700072941 


35 


12084. 


,29442705883 


12084. 


,29442705875 


12084.263 


12083.98 


40 


15785. 


,2626475018 


15785. 


,2626475007 


15784.3 


15785.2626475007 


50 


24667. 


683593322 


24667. 


,683593313 


24667.683593398 


24662.5 



Table 3: The eigenvalues from Example 17. 2^ calculated with the use of the spectral shift technique 
with iV = 40, M = 50000. The spectral shift is given by A^,*^^ = 50n + (2n + 0.5)i. 



mention also that the same problem was considered in [3j where less accurate results were reported: 
Ai = 7.37398502, A2 = 36.3360196, A3 = 85.2925811, A4 = 154.098619 and A5 = 242.705545. 

The SPPS representation allows one to calculate easily the approximate eigenfunctions. The 
eigenfunctions that were obtained for this example are shown in Fig. [2j 

Example 7.3. Consider the following spectral problem [6| Example 2]. 



-y" + + y = Ay, < X < TT 

y(vr,A) = 0. 



(7.5) 



The parameter v was chosen equal to 2. In this example we computed both partial and general 
solutions using the SPPS representations. Since the exact partial solution satisfying the asymptotic 
condition (j2.3p is known and is given by A^/xIl[^) , we compared the approximate solution with 
the exact one. For = 40 and M = 50000 the absolute error was less than 7 • 10"^*^. With the aid 
of Maple software we found that the characteristic equation of problem (|7.5p has the form 

^MA/4,i(vr2) = 0, (7.6) 



where is the Whittaker function, and used equation (j7.6p to compute the exact eigenvalues. 

In Tabled] we present the results obtained by the SPPS method, exact eigenvalues and the results 
from [6j, MATSLISE and SLEIGN2. As in Example 17.11 here again all presented digits in the 
exact eigenvalues coincide with the results delivered by MATSLISE. The performance of the SPPS 
method was considerably better than that of SLEIGN2 and even the 50th eigenvalue was computed 
correctly to 9 decimal places. 

Example 7.4. Consider a particular case of the hydrogen atom equation |6| Example 4]. 

U(^,A) = 0. 
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5 10 15 20 25 30 35 40 45 

Eigenvalue number 



Figure 1: Relative error of the first 50 eigenvalues of the spectral problem for the Boyd equation 
(|73]) calculated with SPPS, SLEIGN2 and MATSLISE. 



Eigenfunctions 

0.3 I 1 1 1 1 1 




X 



Figure 2: The first 10 eigenfunctions of problem (|7.3p calculated with the use of the SPPS method. 



n 




(SPPS) 


^/K (Exact/MATSLISE) 


^/)^ (SLEIGN2) 


([6J) 


1 


2 


4629499739737 


2.4629499739740 


2.46295003 


2.462949030 


2 


3 


2883529299493 


3.2883529299426 


3.28835311 


3.288339398 


3 


4 


1498642187456 


4.1498642187448 


4.14986471 


4.149833151 


4 


5 


0636688237348 


5.0636688237341 


5.0636695 


5.063634795 


5 


6 


0075814581165 


6.0075814581160 


6.0075836 


6.007577378 


7 


7 


9397373768999 


7.9397373768993 


7.939745 




10 


10 


8861250916182 


10.8861250916173 


10.886149 




15 


15 


8426318195682 


15.8426318195682 


15.84275 




20 


20 


8202301908125 


20.8202301908124 


20.82057 




30 


30 


7973502195887 


30.7973502195868 


30.7989 




50 


50 


77867680977 


50.77867680951 


50.789 





Table 4: The eigenvalues from Example 17.31 calculated with the use of the spectral shift technique 
with = 50, M = 50000. The spectral shift is given by given by a[)"^ = lOn + (n + 
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XI (SPPS) (exact/MATSLISE) ^/X^ (SLEIGN2) ^/K ([6j) 



1 


1.97027445061574 


1.97027445061572 


1.9702743 


1.970274439470 


2 


3.00436042551872 


3.00436042551857 


3.0043600 


3.004360435708 


3 


4.01515351791731 


4.01515351791736 


4.015153523 


4.015153641332 


4 


5.0193472218607 


5.0193472218612 


5.0193459 


5.019347630098 


5 


6.0210053515482 


6.0210053515488 


6.0210049 


6.021006315094 


7 


8.0215089715470 


8.0215089715478 


8.0215092 




10 


11.0202653559392 


11.0202653559399 


11.0202663 




15 


16.0176675547289 


16.0176675547294 


16.0176656 




20 


21.0155251794158 


21.0155251794156 


21.0155504 




30 


31.0125189152594 


31.0125189152597 


31.01269 




50 


51.00916429569 


51.00916429551 


51.0097 





Table 5: The values of ^/X^ from Example [731 calculated with the use of the spectral shift technique 
with N = iO, M = 50000. The spectral shift is given by a[,"^ = lOn + (n + 



We have chosen c = 6 and found with the help of Maple the characteristic equation of problem 

5(2i\/A7r) = 0, (7.8) 



8A^A 2v^'2 

where M » 5 is the Whittaker function. In this example we computed both partial and gen- 

eral solutions using the SPPS representations. The procedure for computing the eigenvalues for 
this example is completely analogous to what was described above. In Table [5] we present the 
results obtained by the SPPS method, exact eigenvalues and the results from [6], MATSLISE and 
SLEIGN2. Again all presented digits in the exact eigenvalues coincide with the results delivered 
by MATSLISE. The SPPS method and SLEIGN2 performed similarly to the previous example. 

The next example shows that the situation may change significantly when another value of the 
parameter c is chosen. 

Example 7.5. Consider the same problem as in Example 17.41 but c = —1/4 which corresponds to 
I = —1/2 in (12. 6p . The computation of eigenvalues for problem (17. 7p performed by MATSLISE 
took several hours on Intel 17-3770 microprocessor meanwhile the computation time required by 
SLEIGN2 and SPPS did not change significantly. We present the relative error of the first 50 
eigenvalues on Figure [3l For the SPPS method we used = 40 and M = 1000000. The spectral 
shift was given by a[,"^ = lOn + (n + The accuracy of the results delivered by MATSLISE was 
considerably lower and once again we emphasize the stability of the accuracy of the eigenvalues 
computed by the SPPS method in comparison to SLEIGN2. It is worth mentioning that not only 
for the extreme value / = —1/2 but also for the values close to —1/2 the computation time required 
by the MATSLISE package increases significantly. 

Example 7.6. Consider the following spectral problem [6l Example 3]. 

i-y" + + smx)y = \y, < x < n 

For this problem we were unable to find an exact characteristic equation. We computed both 
particular and general solutions using the SPPS representations. The obtained results are presented 
in Table El The eigenvalues computed by the SPPS method are very close to those delivered by 
MATSLISE. 
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Figure 3: Relative error of the first 50 eigenvalues of problem (|7.7p with c = —1/4 calculated with 
SPPS, SLEIGN2 and MATSLISE. 



n 




(SPPS) 


y/K, (MATSLISE) 




% (SLEIGN2) 


Results from [6] 


1 


1 


69965392162508 


1 


69965392162512 


1 


69965496630163 


1.699674822427 


2 


2 


60438727880122 


2 


60438727880111 


2 


60439344911062 


2.604506077325 


3 


3 


56972957088682 


3 


56972957088910 


3 


56974717663817 


3.570068095387 


4 


4 


55232022604084 


4 


55232022604096 


4 


55235767643976 


4.553053525686 


5 


5 


54189892161891 


5 


54189892161906 


5 


54196736047125 


5.543261224280 


7 


7 


53001773432564 


7 


53001773432606 


7 


53019003621712 




10 


10 


5211087141260 


10 


5211087141255 


10 


5215767716971 




15 


15 


5141539227764 


15 


5141539227760 


15 


5156303086973 




20 


20 


5106568768322 


20 


5106568768319 


20 


5139986174471 




30 


30 


5071385063024 


30 


5071385063018 


30 


5174365365478 




50 


50 


5043027454211 


50 


5043027452760 


50 


5434201591972 





Table 6: The values of ^/X^ from Example 1 7. 61 calculated with the use of the spectral shift technique 
with N = m, M = 50000. The spectral shift is given by a[,"^ = lOn + (n + 
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Figure 4: The graphs of min|;^|=r |<I>7v(A)|, max\x\=r\^W ~ ^nW\ and estimate (I2.16P as an 
ihustration for Rouche's theorem in Example 17.71 for = 50 and = 75. 



7.2 Complex spectrum 

Numerical tests discussed in the previous subsection show that the SPPS method is highly compet- 
itive with the best existing codes on their field of applicability. It delivers stable and reliable results 
even though there remains still plenty of room for improving different computational aspects of the 
developed programs which implement the SPPS method. Moreover, the range of applicability of 
the SPPS method to the difference from the other considered codes includes complex coefficients, 
differential operator on the right-hand side of (j2.6p and complex eigenvalues. Here we present one 
such example. 

Example 7.7. Consider the following spectral problem with the right-hand side of the equation 
involving a derivative 

f-y" + ^y = Ay', 0<x<l 
\y'(l,A) = 0. 

Using Maple we found that the solution of equation (j7.10p satisfying the asymptotic condition (|2.3p 
is given by the expression 



(7.10) 



y{x;X) 



22'+ir(/ + 3/2) 



where / is the modified Bessel function of the first kind. The derivative of this solution has the 
form 

- xJ^y^ [^'^ + (t ) - ^^^'-/^ (t ) 

For a numerical experiment we have chosen / = 1/2, hence the exact characteristic equation is 



y'{x;X) 



0. 



(7.11) 



When a problem admits complex eigenvalues we need to distinguish which roots of the poly- 
nomial <I>Ar(A) correspond to the eigenvalues and which are spurious roots appearing due to the 
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Figure 5: The image above corresponds to the graph of — log |$75(A)|. The semicyhnder marks 
the boundary of the disk |A| =24 illustrating the region of applicability of Rouche's theorem, see 
Example 17.71 The six peaks inside the region represent approximate eigenvalues. The rest of the 
peaks correspond to the roots of the polynomial ^j^^X) appearing due to the truncation procedure 
and do not approximate the true eigenvalues of the problem. 
The image below is the top view of the one above. 



22 



n 



Xn (SPPS) 



A„ (Exact) 



1 4.47123493365+ 6.76481747492i 

2 5.63553225528+ 13.37799928406i 

3 6.35749327967+ 19.82515033089i 

4 6.88515096019+26.20887598267? 

5 7.30184486323+ 32.5608828157H 
10 8.62739882797+ 64. 14303168943i 
20 9.98333956705+ 127.0816376255? 
30 10.7844002548+189.9555609957? 
50 11.7983559318+315.6569255418? 
75 12. 6055406455+472. 7574366522i 
100 13. 1790674123+ 629.8482784849i 



4.47123493371+6.76481747480? 
5.63553225515+ 13.37799928396? 
6.35749327947+ 19.82515033081? 
6.88515095992+26.20887598266? 
7.30184486294+32.56088281579? 
8.62739882786+64.14303168978? 
9.98333956726+ 127.0816376257? 
10.7844002552+ 189.9555609955? 
11.7983559297+315.6569255437? 
12.6055406452+472.7574366509? 
13. 1790674160+ 629.8482784850i 



Table 7: The values of A„ from Example 17. 71 calculated with the use of the spectral shift technique, 
M = 200000 and N = 50. 



truncation procedure. Contrary to the case of a purely real spectrum we cannot simply discard 
all roots whose imaginary part is greater than some e. Instead, Rouche's theorem and estimate 
()2.16p suggest that the roots closest to the origin (or to the current centre Aq when the spectral 
shift is used) cannot be the spurious roots. We give an illustration of application of this theorem. 
According to Rouche's theorem we need to find such values of the radius r that 

min |«>Ar(A)| > max|^>(A) - $7v(A)|, (7.12) 

|A|=r \M=i" 

establishing that the numbers of zeros of the functions ^n{X) and ^*(A) coincide inside the disk 
|A| < r. We calculate min|;^|=r |*^*Af(A)| using the SPPS representation. To estimate the difference 
|<1>(A) — $7v(A)| we use (j2.16p . According to Lemma [2m for problem (j7.10p one can take C = 3/2. 
We present the obtained results on Figure H] for A'^ = 50 and N = 75 where it is compared to 
max|;^l=^ |^('^) ~ *^*Af(A)| calculated with the aid of the exact characteristic function ()7.1ip . As can 
be seen from the graphs, when we use the exact characteristic function for the estimation of the 
error the largest values of r for which inequality (j7.12p holds are rso ~ 17 and rys ~ 25. The 
estimate (I2.16P delivers rougher estimates, r^g ~ 7 and r^g ~ 12. On Figure [5] we show the graph of 
— log |<1>75(A)| together with the boundary of the disk |A| = 24. The peaks on the graph correspond 
to the roots of the polynomial $75. Only six of them are located in the disk and hence should be 
considered as approximations to the eigenvalues. All other roots are located outside the disk and 
hence should be discarded. 

To calculate the higher order eigenvalues, we applied the spectral shift technique described in 
Section [5l The following strategy for choosing the values of the spectral shift was implemented. 

(n) 

Let A; be the spectral shift on the n-th step and Ai, . . . , A„ the already found eigenvalues. Based 
on the representations (|5.7p and (|5.8p we calculate the roots of the polynomial $^^(A) and reorder 

(n) 

them with respect to the distance from A^ . From these ordered roots we choose the closest 

(n) 

to the point A* and sufficiently distant from Ai,...,An. We denote this root as A,i_|_i and set 
^(n+i) _ _|_ where AA is a fixed displacement. Performing this procedure with a1*^^ = 0, 
AA = — ?, = 50 and M = 200000 we obtained the list of eigenvalues presented in Tabled The 
eigenvalues are ordered according to their distance to the origin. As can be seen from the table the 
approximate eigenvalues are calculated with a remarkable accuracy. 
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